Mon. Not. R. Astron. Soc. 000, [[ju] (2002) Printed 1 February 2008 (MN M^X style file v2.2) 



A general, three-dimensional fluid dynamics code for stars 
in binary systems 



(N 
O 
O 
(N 

u 

< 
in 

(N 

> 
oo 

o 



6 



Martin E. Beer* and Philipp Podsiadlowski 

University of Oxford, Wilkinson Building, Keble Road, Oxford, OX1 3RH, England 



Accepted 2002 April 22. 



ABSTRACT 

We describe the theory and implementation of a three-dimensional fluid dynamics 
code which we have developed for calculating the surface geometry and circulation 
currents in the secondaries of interacting binary systems. The main method is based 
on an Eulerian-Lagrangian scheme to solve the advective and force terms in Euler's 
equation. Surface normalised spherical polar coordinates are used to allow the accurate 
modelling of the surface of the star, as is necessary when free surfaces and irradiation 
effects arc to be considered. The potential and its gradient are expressed as sums of 
Legendre polynomials, which allows a very efficient solution of Poisson's equation. The 
basic solution scheme, based on operator splitting, is outlined, and standard numerical 
tests are presented. 
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1 INTRODUCTION 



Irradiation of the secondaries in X-ray binaries can dramat- 
ically change their appearance and their internal structure. 
The irradiation pressure force can lead to significant distor- 
tions of the surface (Phillips & Podsiadlowski 2002), while 
irradiation-driven circulation currents can transport signif- 
icant amounts of energy to the unirradiated side. There is 
ample observational evidence for the existence of such cir- 
culation currents: e.g. in HZ Her/Her X-l (Kippenhahn & 
Thomas 1979; Schandl, Meyer-Hofmeister & Meyer 1997), 
in cataclysmic variables (Davey & Smith 1992) , in Nova Sco 
during outburst (Shahbaz et al. 2000) and in Cyg X-2 (J. 
Casares & P. A. Charles 1999; private communication), there 
is clear evidence that a substantial amount of X-ray heated 
material can move beyond the X-ray horizon. Moreover, per- 
sistent residuals in the observed radial-velocity curves in X- 
ray binaries, e.g. in Nova Sco during outburst (Orosz & Bai- 
lyn 1997) and Vela X-l (Barziv et al. 2001), may provide 
direct evidence for circulation. 

Modelling of the irradiation-induced circulation in bi- 
naries is difficult due the three-dimensional nature of the 
problem. It requires the simultaneous solution of the shape 
of the irradiated star and the circulation and a proper treat- 
ment of the surface boundary conditions. 

Attempts to model circulation in irradiated secondaries 
have been made by various authors in the past, initially 
using perturbative methods either in planar geometry (Kip- 
penhahn & Thomas 1979; Kirbiyik & Smith 1976) or spher- 
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ical geometry (Kirbiyik 1982), and more recently non- 
perturbatively using smooth-particle hydrodynamics (Mar- 
tin & Davey 1995). None of these investigations to-date, 
however, are self-consistent, neither allowing for changes in 
the surface geometry nor including radiative surface stresses. 

On the other hand, circulation currents have been ex- 
tensively modelled in multi-dimensions in the context of 
modelling the circulation in the ocean and tidal flows on 
Earth, where over the years efficient methods have been de- 
veloped to treat circulation realistically and in a numeri- 
cally efficient way. A commonly used method is based on 
an Eulerian-Lagrangian scheme (a description of which may 
be found in Lu & Wai 1998). This method is similar to the 
upwind method but more physically sound in the physical 
treatment of advective terms and has been shown to be un- 
conditionally stable (Casulli 1990; Casulli & Cheng 1992). 

The purpose of this paper is to present the philosophy 
and the details of a fairly general three-dimensional fluid dy- 
namics code which we have specially developed for treating 
the secondaries in interacting binaries, in particular under 
the influence of external irradiation. The main method is an 
application of the Eulerian-Lagrangian method by Lu & Wai 
(1998) which we have modified for our application, draw- 
ing also on the results of related work by Uryfl & Eriguchi 
(1995, 1996, 1998), Miiller & Steinmetz (1995) and using 
methods developed in the context of geophysical fluids (see 
e.g. Pedlosky 1987). At present our code is still somewhat 
simplified since we do not include the thermodynamic equa- 
tion, necessitating the use of a polytropic equation of state. 
In Appendix B, we describe how we plan to generalize the 
code in the future. In a subsequent paper, we will apply the 
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Figure 1. Schematic diagram defining the coordinate systems (r, 6, <j>) and (x, J/, 2) relative to the system geometry and other binary 
parameters (A: orbital separation, My. mass of the primary, f2: angular velocity of the secondary). 



code first to rotation in the standard Roche problems, and 
then to study irradiation-driven flows in X-ray binaries, con- 
sidering both the effects of heating and external radiation 
pressure, and of radiative surface stresses. 

The outline of the paper is as follows. In Section | we de- 
scribe the transformations and the usefulness of surface fit- 
ting coordinates, Section ^ describes the basic equations and 
the dimensionless variables used. The theory and implemen- 
tation of the calculation of the gradient of the gravitational 
potential is given in Section ||, including estimates of its 
accuracy. The general solution method is described in Sec- 
tion p. Finally, Sections ||and§ present standard numerical 
tests of the code, advection tests and Maclaurin spheroids. 



2 COORDINATE SYSTEM 

Figure |l| defines the adopted coordinate system, standard 
spherical polar coordinates centered on the center of mass 
of the secondary, where the directions of the axes are given 
by the unit vectors (r, 8, <p). 



2.1 Surface normalised coordinates 

For our main applications it is important to model the stel- 
lar surface accurately and to allow it to adjust freely to 
satisfy whatever surface boundary conditions are applied. 
If the outer edge of the coordinate grid did not coincide 
with the surface, it would be non-trivial to calculate surface 
stresses and derivatives along the surface accurately. These 
difficulties can be avoided by using a grid whose boundary 
is defined by the stellar surface. To achieve this, we follow 
Uryu & Eriguchi (1996) and transform our basic equations 
(see Section []) using spherical surface fitting coordinates, 
i.e. 



(1) 



where R(0, <f>) is the radius of the star in the direction (9, </>). 
With this definition, r* is restricted to the range 

< r* < 1 . (2) 



This means that the grid has to adapt continually as the 
stellar surface changes and that derivatives are transformed 
according to: 
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(4) 
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(6) 



The surface derivatives J^- and J^pr are calculated using the 
potential as described in Appendix C3. 



2.2 Coordinate grid 

To define the coordinate grid, we split the star into equally 
spaced intervals in r, cos 9 and <j> with a total of N T , N$ , N$ 
elements. The equal spacing in cos# ensures that all grid 
cells have the same volume. The centers of the grid cells in 
the 9 and <f> directions are defined by 



cos 9\ 



U - 0-5) 



2n(k - 



(7) 



where we restricted the 9 range assuming even symmetry 
with respect to the xy-plane, and the surface normalised 
radial coordinate becomes 



(^-0.5)^(^,00 



0.5 



Nr 



(8) 



Velocities are calculated on cell boundaries whilst pressures 
and densities are evaluated at the centres of cells. Figure ^ 
shows a cross-section of the grid indicating where various 
quantities are evaluated. 

This grid spacing acts so that the entire grid adapts 
when the surface changes. We adapt this so that if desired 
the central grid region stays fixed at its initial position when 
the surface adjusts. This enables us to fix the innermost 
regions of the grid in situations when we are only interested 
in the outermost layers. For this purpose we define N c as 



A three-dimensional fluid dynamics code for stars in binaries 3 









*Pi-Uk 


„ K, j>k 

3 — • C 


' i,j,k 

3 — • 

^'Pi + l,j,k 




"a 
. ij+l.k 

W\ J+ l.k 





Figure 2. Cross-section of the adopted coordinate grid in the r 
and 8 directions, indicating where quantities (velocities, u, pres- 
sure, p, and density, p, are evaluated. 



the number of radial layers from the center whose positions 
are kept fixed. 



(i-0.5)flin(9j,4> k ) 

JVr 

,,. _ , JV c fl in (9j,0 k ) + (»-0.5-iVc) x 



(Ni—N ) 
jV c fl in (9j,0 k ) 
AT, 



if i < iV c 



otherwise 



(9) 



If desired it is possible to redefine the radial spacing so that 
it is no longer linear but concentrated in a particular region, 
e.g. near the surface. The number of grid points, however, is 
kept fixed during a calculation. Fixing the innermost region 
of the grid also requires an alteration to equations (J3| — , 
since the term proportional to the surface curvature has to 
be replaced by a term representing the curvature of the grid 
at the required radius. 
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if i < N c 



otherwise 



(10) 



3 BASIC EQUATIONS 

The equation of motion relating velocity (it), pressure (p), 
density (p), frictional force (^F) and potential ($') in a frame 
rotating about the centre of mass of the secondary with 
angular velocity fi is 

^ + 2Qxti = -ivp- V$' +T . 
at p 

where 
du du 
~dt ~ ~dt 

and the potential term is given by 

GMxr sin 6» cos (fixr) 2 



+ (u-V)u 



(11) 



(12) 



$' = $ + 



A 2 



GMi 

r-2 



(13) 



where the second term represents a coordinate transforma- 
tion from the centre of rotation to the centre of the sec- 
ondary. Here Mi is the primary mass, G the gravitational 
constant, A the binary separation, r2 is the distance from 
the primary to the point (r, 8, <j>), and the potential of the 
secondary $ is given by Poisson's equation 



= 4n Gp . 
The continuity equation, 
dp 



dt 



+ pV-u = 



(14) 



(15) 



can be written in integral form, using the divergence theo- 
rem, as 

d 



pu-n dS = 



(16) 



where n is a unit surface normal vector. These equations 
need to be solved along with an equation of state, assumed 
to be a polytropic equation of state at present which relates 
pressure and density according to 



Kp 



(17) 



where K and n are the polytropic constant and polytropic 
index, respectively. 

3.1 Boundary conditions 

To completely specify the mathematical problem, we also 
need to specify a set of boundary conditions at the centre 
and the surface. The central boundary conditions are un- 
problematic and are given by 

du T 
dr 

The simplest boundary condition for the surface are the 
zero-pressure condition and the conservation of the mass 
of the secondary, 



u = 



= 



(18) 







M2 = constant , 



(19) 



(20) 



or, more realistically, a pressure condition that assumes an 
atmosphere in hydrostatic equilibrium, i.e. 

2<7 

P= it , 
6k, 

where g is the effective surface gravity, and k the photo- 
spheric value of the opacity (see e.g. Kippenhahn & Weigert 
1994). 

In the case of irradiated stars, a special treatment of 
the surface layer is required since the irradiation flux is de- 
posited in a thin (turbulent) surface layer (of order an atmo- 
spheric scale height), which will generally be smaller than 
our grid size. A separate atmosphere calculation, which in- 
cludes the effects of heating and irradiation pressure, then 
determines the pressure across the surface of the star (anal- 
ogously to the case of normal stars; also see Tout et al. 
1989). In addition, the variation of the radiation pressure 
force causes a surface stress which drives horizontal motion 
perpendicular to the stress and vertical motion in a thin tur- 
bulent boundary layer (an 'Ekman' layer), a process known 
as 'Ekman' pumping. This produces a vertical velocity com- 
ponent in regions where the surface stress varies. This pro- 
cess is entirely analogous to the wind-driven circulation in 
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oceanic circulation systems (see chapter 5 of Pedlosky 1987) 
and can be treated analogously. 

3.2 Surface adjustment 

At the end of each time step, the surface is adjusted using 
the current values of the velocity at the free surface. The 
surface normal is given by 

. . 18R- 1 8R Z 

n = r ^~~^ e ^77J^ ■ ( 21 ) 

r 08 rsmff a<p 

The dot product of the surface velocity with the surface 
normal gives the velocity component normal to the surface. 
Equating this to the component of the surface adjustment 
in the radial direction (Sr), which is normal to the surface, 
yields 



Sr 



u-n 
n-f 



(22) 



This is added to the current value of r at the end of each 
time step. In order to ensure mass conservation, the density 
of each 8, <f> element is scaled by the ratio of the old and new 
radii cubed 

H^) 3 "- (23) 

This ensures that the density decreases as the radius in- 
creases, conserving mass in the process. If mass is not con- 
served instabilities can develop at the surface of the star. For 
example, consider the case of a contracting star. As mass 
leaves the inner radial boundary of the surface layer, no 
mass loss/gain occurs through the upper boundary and so 
the density in the surface layer drops to zero. The central 
density is calculated by extrapolating the densities in the 
central regions of the star. It is only used as a boundary 
condition in the potential calculation (Section ^]). The ex- 
trapolation occurs by extrapolating in the radial direction 
for each 8 and cf> direction before taking the mean of these 
values. 



3.3 Viscosity 

Artificial viscosity is often added in numerical simulations 
to smooth out the flow and to broaden shock fronts. Stone & 
Norman (1992) in their code ZEUS-2D use an artificial pres- 
sure similar to that of von Neumann & Richtmyer (1950): 



1 = 



C 2 p(Av) 2 




if (An) < 
otherwise , 



(24) 



where Av is the change in velocity across a cell and C = 
I /Ax. C measures the number of zones over which the ar- 
tificial viscosity spreads a shock and is typically chosen to 
be C ~ 3. This has been generalised to multidimensions 
by Tscharnuter & Winkler (1979) using an artificial viscous 
pressure tensor (Q) 



Q = 



l 2 pV-u [Vu - |(V-«)e] 







if V-u < 

otherwise , 



(25) 



where e is the unity tensor, Vu is the symmetrized tensor 
of the velocity field and I is of the order of the local width 
of the grid. Following Stone & Norman (1992), we neglect 
the off-diagonal (shear) components of the artificial viscous 



pressure tensor. The (nonzero) diagonal elements of the grid 
are 



Qn = I pV-tt 



Q22 = I pV-u 



Q33 = I pV-u 



du r 1 

■oF-3 (v ' u) 



1 dug U r 1 

r 06 r 3 



1 Bus u r Ug 

— -JUT H 1 cot ' 

r sm 8 o(p r r 



(26) 
(27) 



V-u) 



(28) 



The frictional force in the momentum equation ([ll]) is 
_ 1 



-V-Q 



where 

- , 1 <9Qii 
r.T = — 



(29) 



e.T 



3Qn 

pr 

(Q11 + 2Q 22 ) 



p Or 

1 dQ 2 2 
pr 88 pr 

1 d(Qn + Q22) 



cot 8 , 



(30) 

pr sm U Ocp 

and where we have used the property Tr(Q) = to eliminate 
Q33. This source of viscosity may be switched on or off in 
the code. 

Stone & Norman (1992) also consider an artificial linear 
viscous pressure, 



q = CipC^Av , 



(31) 



to damp instabilities in stagnant regions of the flow, where 
Ci is a constant of order unity and C a is the adiabatic sound 
speed. This is calculated separately for each direction and 
then added, e.g. in the r direction, as 



p or 



(32) 



In our code we include the linear viscous pressure which may 
be switched on or off. 



3.4 Dimensionless variables 

In the code we use dimensionless variables defined as: 



R 



R_ 

R B 



P - 



P_ 

Pc 



$ = 



M = 
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GpcRs 2 ' 


RsVUpl 


M 




Lx 




GR 2 c ' 


T 


V - 


V 


Gp c R s ' 


RJVGpl 



Q = 




R = Kp c R s 



(33) 



where _R S is the solar radius, p c the initial central density 
and G the gravitational constant. This transforms equa- 
tions (O), (M) & (0) to: 
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(34) 
(35) 

(36) 



^ + 2ft x u = -~Vp - V$' + T , 

dt p 

/ pdV + / pu-ndS = , 

9* ,/dv 

P GR^pS 1 -^ ' 

4 POTENTIAL CALCULATION 
4.1 Theory 

Miiller & Steinmetz (1995) developed an efficient algorithm 
for solving Poisson's equation which utilizes spherical co- 
ordinates and an expansion into spherical harmonics. This 
results in an algorithm which has a computational cost pro- 
portional to (L + l) 2 NrNoNj, where L is the highest order 
harmonic considered. The general solution of Poisson's equa- 
tion in spherical harmonics, Y lm (8,(f>) (Morse & Feshbach 
1953) can be written as 

oc I 

<E.(r,e,0) = -G^^ T ]T Y lm (8,^)x 



p(8' ,</)' ,r') sin 8' Pi m (cos 8') is an odd function and its inte- 
gral from 8 — to ir in equations ( p^ ) & ( |39| ) is zero. Hence 
terms odd in I + m do not contribute to the sum, and we 
only need to include values of m from to I and double the 
contribution of the positive m term. 

In our code we deal with $' rather than 



*'(r,M) =*(r,M)- 

(S>sin#) 2 GAfirsin6»cos0 GM\ 



A 2 



A 2 



(A 2 + r 2 — 2Ar sin 8 cos <f>) 2 



GM 1 r cos 6* cos 



(rf2) sin cos #+ 



GM\Ar cos # cos ( 



,4 2 



(A 2 + r 2 - 2Ar sin 8 cos 0) 5 



(44) 



This then yields: 

d& 9$ a 2 ■ 2 a^ 
-r— = — r\l sm 

ar ar 

GM\ sin # cos GM\ (r — A sin 8 cos 0) 

+ ~ — : : : — > (45) 



(46) 



9$' d$> 



1 I T^lm / \ 

+ r D (r) 



(37) 



where 



C lm {r) 



D lm (r) 



/ dftV 

J4tt 



m *{8'A') I dr'r' l+2 p(r',8'^') , (38) 



dQ'Y lm *{8',(f>') I dr'r n ' l p{r',e',cj)') , (39) 



dQ,' = sin d'dd'dtj)' 



(40) 



By differentiating this equation we may obtain analytical 
formulae for the gradient of the potential in spherical coor- 
dinates (as derived in Appendix CI): 



d$(r,6>, 
dr 



+ 



1=0 m=—l 

r i +2 -° (r)+lr D (r) 



(41) 



d$(r,fl, 
88 



,|m+l| 



(cos 8) mcost 



P ; (cos 
d$(r,8, 



sin # 



■J+iC (r) + rD (r) 



(42) 



47T 

1=0 m=-l 



[^T^(r) 



+ r'-D !m (r) 



(43) 



For odd i + m, P; m (cos#') is an odd function in 8'. 
Since the density is an even function w.r.t. 8' (because 
of the assumed symmetry with respect to the xy-plane) 



GM\r sin 8 sin < 
A 2 



GM^Ar sva8 sin <f> 
(A 2 + r 2 - 2^4r sin 8 cos cj>) i 



(47) 



4.2 Accuracy 

In Appendix |C2| , we describe the implementation of this po- 
tential calculation in the code. To test the accuracy of the 
potential calculation, we compared it to three simple cases in 
which analytical solutions exist. We also compared it to the 
ZEUS-2D code (Stone & Norman 1992) for which the same 
tests have been performed. For the purposes of comparison, 
we chose their test cases in spherical polar coordinates in 
which they assumed equatorial symmetry as in our calcula- 
tions. They do not list the number of grid points they used. 
Here we used 50 3 grid points. Table [l] shows the comparison 
in the potential calculation for three cases of a homogeneous 
sphere, a centrally condensed sphere and a homogeneous el- 
lipsoid. 

Table [l] shows that our potential calculation is more 
accurate than that of ZEUS-2D regardless of the density dis- 
tribution. As our calculation is based on a numerical repre- 
sentaion of the analytical value it is extremely accurate for 
the simplified case of a homogeneous sphere. The centrally 
condensed sphere shows the accuracy of the radial integra- 
tion which is resolution limited. The error in the calculation 
for the homogeneous ellipsoid is smaller than that of the 
ZEUS-2D calculation. Indeed, it is only large at the surface 
- throughout the remainder of the star it remains of sim- 
ilar order as in the calculation for the centrally condensed 
sphere. We have also calculated the error in the 8 derivative 
of the potential. We find that this is of similar order to the 
error in the r direction. 
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6 M. E. Beer and Ph. Podsiadlowski 



Object 


Max. error in ^ 

or 


Max. error in V$ in ZEUS-2D 


Homogeneous Sphere 


1.02xl0~ 8 % 


i.iixicr 2 % 


Centrally Condensed Sphere 


6.35xl(r 2 % 


9.98xl0- 2 % 


Homogeneous Ellipsoid 


1.53% 


1.86% 



Table 1. A comparison of the potential calculation in our code, performed on a three-dimensional grid, to the two-dimensional calculations 
using the ZEUS-2D code for three simple cases with analytic solutions. 



5.1 Operator splitting of the equation of motion 

The equation of motion is solved in two parts using the 
method of operator splitting or fractional step. With this 
technique the equation of motion is split into two parts. One 
representing the advective terms, and another representing 
the force terms. The advantage of this method is that dif- 
ferent numerical techniques may be used to solve the two 
equations representing physically different processes. 

We split equation ( |ll| ) into two parts, one containing 
the advective terms, the other the force terms. Once the first 
part is solved, its solution is used in the second part to find 
the solution corresponding to the full time step. Represent- 
ing the time derivative of these equations in finite difference 
form yields 

a t 

— ^+2ftxu a = 0, (48) 

ot 

^- = -iv P -V<£-' + ^, (49) 

ot p 

where we have followed the notation of Lu & Wai (1998), 
and the subscript p refers to the position of the element of 
interest in the previous time step. Each part of the equation 
is solved for a full time step {St), and we have used the 
notation u a to indicate the solution for the velocity once 
the first equation, containing the advective terms, has been 
solved. 

By solving the equations using the velocity of the fluid 
element from the previous time step, the non-linear term 
(u-V)u is eliminated from the equations. This means that 
advective terms in the equations are solved in a Lagrangian 
frame (for fixed mass elements) rather than an Eulerian 
frame (with fixed positions) where the remainder of the 
terms are solved for. Hence this method is known as an 
Eulerian-Lagrangian method. Figure ^| illustrates how u* is 
calculated. The figure shows a cut away element of the fluid. 
Using the velocity we calculate the position of the fluid 
element at a time step St before the present one. Interpolat- 
ing the velocity grid yields the velocity tip at this position. In 
actuality, we split the time step into a number of substeps 
enabling an accurate calculation of the streamlines of the 
fluid. The substep method consists of finding the velocity at 
a point a substep away and using this velocity to find the 
velocity in the next substep. This way curved trajectories 
may be easily followed. 

5.2 Solution of the continuity equation on an 
adaptive grid 

The grid used in the code is adaptive. At the end of each 
time step, the grid is rescaled so that its boundary repre- 
sents the surface of the star. Thus, the mass within each cell 




Figure 3. Diagram showing the path interpolation used for the 
Eulerian-Lagrangian method. 

will change at the end of each time step. To compensate for 
this in solving the integral form of the continuity equation 
(equation jl6|), the surface velocity is calculated at each it- 
eration of the continuity equation and the radial velocities 
used in the calculation are relative to this 

u! r (e,<f,) = ur(e,<fi)- 1 - «r(M)l , (so) 

iVr — iV c la 

where the prime indicates the velocity used in the calculation 
and the subscript s a quantity evaluated at the surface. Dorfi 
shows in LeVeque et al. (1998) [p. 279] that time centering of 
the equations results in second order accuracy in St. Hence 
the variables are evaluated at half-odd-integer time steps i.e. 
when calculating the mass flow through a cell boundary we 
use u t+ ~i and p t+ ~i . 

5.3 Calculation of the pressure gradient at the 
surface 

The radial velocity needs to be calculated at the surface 
which would require a ghost point for pressure half a ra- 
dial grid zone past the surface. Instead we use the surface 
pressure which is one of the boundary conditions. Thus we 
evaluate the pressure gradient one quarter of a zone below 
the surface. However, this still yields an incorrect value for 
the surface pressure gradient, so that in equilibrium it will 
not balance the potential gradient causing a spurious sur- 
face velocity. To correct for this, we calculate the potential 
gradient a quarter of a zone below the surface and find that 
this accurately balances the pressure gradient. 
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Setup initial solution or read in model from file 



Use the Eulerian-Lagrangian method to find u a 



Calculate the new value of u using the current 
values of pressure and the potential gradient 



Calculate new values for the density 



Check for convergence of u and p 



Surface adjustment 



Check for convergence of solution 



Write output to file 



Figure 4. A flow chart diagram showing the iteration procedure 
used in the code 

5.4 Iteration procedure 

Figure |^ shows the iteration procedure used in the code. 
After initialization the main loop of the code is entered. The 
Eulerian-Lagrangian method is used to find the solution of 
the advective terms. A sub loop is then entered in which the 
velocities and density are solved. Using the current values 
of pressure and density, the new values for the velocities are 
calculated. With these velocities the new densities can be 
determined. Once these have converged, the end of a time 
step has been reached and the surface is adjusted to fit the 
new shape of the star. Overall convergence is then tested 
until a solution is found which is written to an output file. 
It is not guaranteed that the iteration procedure will always 
converge. 



6 ADVECTION TESTS 

To test the advection in our code, we use two simple tests. 
Both are in spherical geometry and have also been carried 
out by Stone & Norman (1992) for their code ZEUS-2D to 
which we compare it. 

6.1 Expansion of a homogeneous sphere in a 
velocity field 

The first test is a "relaxation" problem for a pressure-free 
and gravity-free homogeneous gas with a velocity field pro- 
portional to the radius (u r = vor). The density will decay 
exponentially with time as 

p(t) = P0 e- 3 ^ . (51) 

At vot = 4, the density should have decayed by almost six 
orders of magnitude to p = 6.144X 10 _6 po- In our calcula- 
tion we use wo = 1, po = 1 and three time steps of 10~ 2 , 1CP 3 
and 1CT 4 . The results are shown in Table 0. 

At t = 4 our code is in better agreement than ZEUS- 
2D (5.60 X 10~ 6 , 8.8%) for all the time steps considered, the 



Time step Density at t = 4 Error in the density at t = 4 

10" 2 6.522 XlO" 6 6.14% 

10" 3 6.181 XlO" 6 0.60% 

10" 4 6.148X10" 6 0.06% 



Table 2. A comparison of the calculated density to the analytical 
solution for the expansion of a homogeneous sphere in a velocity 
field. 



Time step 


Max. error in density at t = 0.535 


lO" 3 


5.30% 


io- 4 


0.50% 



Table 3. A comparison of density calculated in the code to the 
analytical solution for the collapse of a pressure-free sphere. 

error in the density being directly proportional to the time 
step considered. Stone & Norman (1992), however, provide 
no information on the time step they use, so this cannot be 
compared further. 



6.2 Pressure-free collapse of a sphere under 
gravity 

The second advection test we use is the collapse of a homo- 
geneous, pressure-free sphere under gravity. Hunter (1962) 
showed that 

— = cos 2 (3 , — = cos' 6 13 , (52) 
ro Po 

where 

p+^=t ] m^, (53) 

and ro and po are the initial values of the radius and the 
density, respectively. For ro = po = G = 1, the free-fall time 
(the time at which the sphere has collapsed to a point) is 
0.543. The test was run at two different time steps (10 -3 
and 10 -4 ), and the density, radius and the mass were com- 
pared to the analytical solution and the results obtained 
with ZEUS-2D. 

Table ^ shows the errors in the calculated density in 
comparison to the analytical solution. In all cases, the 
density profile remained flat and the mass constant. This 
demonstrates the ability of the code to adjust the surface of 
the star correctly. Figure |B| shows a comparison at time t = 
0.535 with the analytical solution for a time step of 10 -4 . 
Note how the grid is only defined within the sphere allowing 
an accurate representation of the problem with no spurious 
points exterior to the surface. In the calculation with ZEUS- 
2D, a number of grid points end up outside the surface - a 
direct consequence of their use of ghost zones beyond the 
boundary. 
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a 
a 




Figure 5. Test simulation for the pressure-free collapse of a ho- 
mogeneous sphere under gravity. The comparison shows the den- 
sity as a function of radius at dimensionless time t = 0.535 for 
a time step of 10 — 4 . The solid curves represent the analytical 
solution and the open circles our model results. 



7.1 The compressible Maclaurin sequence 

Lai, Rasio & Shapiro (1993, hereafter LRS) calculated a se- 
quence of compressible Maclaurin spheroids based on an en- 
ergy variational method. LRS calculate universal dimension- 
less quantities that are functions of the eccentricity of the 
spheroid only. The eccentricity is defined as 



r>2 

-""pole 
equator 



(54) 



where R po \ e and i? cqU ator are the radii of the pole and equa- 
tor respectively. Transforming into the dimensionless units 
used in our code, their equation (3.27) becomes 



9/ 



K n {l — n/5) 



p c 



where 



P = 



M 



is the mean density and 

- 2 — 

R (^polc-^cquator ) ^ > 



(55) 



(56) 



is the mean value of the radius of the star. 

For a polytropic system, as the polytropic index, n, in- 
creases, so does the radius of the polytrope (which is defined 
as the position where the density drops to zero). For non- 
zero n only one analytical solution exists which has a finite 
radius. This is the case n — 1 and has the solution 

sinf 

where 



(57) 



AttG 



(n+l)K' 



(58) 



This is one of the values for the polytropic index considered 
in both the calculations of LRS and Uryu & Eriguchi (1998). 
Other values for the polytropic index of interest are n = 
1.5, 3. These polytropic indices are reasonable representa- 
tions for the internal structure of convective and radiative 
stars, respectively. Table ^ gives the ratio of mean density 



n 


P/pc 


1 


3.03963 XlO" 1 


1.5 


1.66925 XlO" 1 


3 


1.84561 X10~ 2 


5 


0.0 



Table 4. Values of £ at the surface of a polytrope and the ratio 
of mean density to central density for various polytropic indices 
n. 



o 
■+J 

c 

O 
O 
H 




Figure 6. Angular velocity squared against eccentricity of an 
ellipsoid. The dashed curve represents the results of LRS and the 
solid curve our results for comparison. 



to central density (p/p c ) for various polytropic indices (Kip- 
penhahn & Weigert 1994). 

In our tests we considered the case n = 1. This is the 
only analytical solution which has a finite radius making 
it a convenient starting approximation. Our code can also 
run, however, with approximate solutions appropriate for 
n = 1.5 & 3 using the method of Liu (1996) to calculate the 
density distributions. Starting from this analytic solution 
we have calculated a compressible Maclaurin sequence for 
comparison, where we use 50 grid zones in the r direction 
and 48 grid zones in the 9 direction. We increased the value 
of Q in units of 0.1 and used a time step of 5t — 0.001 to 
find a converged solution. We also included a linear artificial 
viscous pressure to damp out oscillations. 

Figure ^| represents a comparison to the calculations of 
LRS. It shows a plot of the dimensionless angular velocity 
squared against eccentricity. The figure shows that the code 
accurately calculates the eccentricity for each value of fi. 



8 FUTURE APPLICATIONS 

Our first application of the code will be the standard Roche 
problems. These have recently been considered by Uryu & 
Eriguchi (1998), who assumed irrotational flows. We will ap- 
ply our code to the same case for comparison and then relax 
the constraint of irrotation, using realistic surface bound- 
ary conditions. The second and main application will be the 
case of irradiation by X-rays in interacting binaries, where 
we will study the effects of surface heating, irradiation pres- 
sure and surface stresses to obtain self-consistent solutions 
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of the geometry and the circulation system in irradiated sys- 
tems and compare these results to observed systems where 
irradiation effects have been identified (e.g. HZ Her/Her X- 
f , Cyg X-2, Nova Sco, Vela X-l, Sco X-f). At a later stage 
we will extend the code to more realistic stellar models as 
outlined in Appendix H. 
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APPENDIX A: SYSTEM VARIABLES 

Following are a list of variables used in this paper and a description of what they represent: 



d_ 

,11. 



Ill 



Lagrangian derivative 
Eulerian derivative 



p Pressure 

p Density 

T Frictional force 

u Velocity in the rotating frame 

u a Velocity aft er the solution of the advective terms in the equation of motion 

(section |5.l| ) 

fl Angular velocity of rotation 

$ Potential 

G Gravitational constant 

g Surface gravity 

r, 9, cf> Spherical Polar coordinates 

A Separation of the binary 

T2 Distance from the primary to the point (r, 6, (f>) 

M\ Mass of the primary 

M2 Mass of the secondary 

n Unit surface normal vector 

r, 8, </> Unit vectors in the spherical coordinate system 

K, n Constants in the polytropic equation of state (equat ion p. 71 ) 

Superscript denoting dimensionless variable (section p. 4| ) 
N r , Ng, Nff, Number of r, 9, cf> elements 

N c Number of innermost radial zones whose position is kept fixed 

i,j,k Subscripts used to denote r, 9, elements 

t Local optical depth in the secondary 

k Local value of the opacity 

L x X-ray luminosity of the compact object 

q Artificial viscous pressure 

Q Artificial viscous pressure tensor 

Ci Constant of order unity used in the calculation of artificial viscous pressure 

C a The adiabatic speed of sound 

£, ui Variables in the Lane-Emden equation for a polytrope (section 

V Volume of the secondary 

7 Moment of inertia of the secondary 

J Angular momentum of the system 

T Kinetic energy of the secondary 

W Potential energy of the secondary 

Superscript denoting universal dimensionless variable with no dependence on 

polytropic index (section 

R Mean radius of the secondary 

_R P oie Polar radius of an ellipsoid 

^equator Equatorial radius of an ellipsoid 

Ro Radius of a spherical polytrope of mass equal to the secondary 



APPENDIX B: GENERALISATION TO A NON-POLYTROPE 

In this appendix we describe how the computer code may be extended to more realistic stellar models. The most important 
change is to use a realistic (non-polytropic) equation of state and the addition of an energy equation, the thermodynamic 
equation (Tassoul 1978): 

p^+pV-it = *„ + /9e nuc - V-(F + F r ) , (Bl) 

where U is the total internal energy per unit mass, the heat generation by viscous friction, e nuc the rate of energy released 
by thermonuclear reactions per unit mass, F the heat flux, 

F = - X VT , (B2) 



and F T is the radiative flux, 
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Fr = _*acJ_ VT (B3) 

where a is the radiation constant and T the temperature. The thermodynamic equation may also be written as 

pT ^- =*„ + pe nuc - V-(F + F r ) , (B4) 
at 

where S is the entropy per unit mass. Landau & Lifshitz (1987) show that 



dt ( dt) p dt ' vs ( at ) p VT 

Therefore, 



(B5) 



T ^ = Cp dF' (B6) 
where c p is the specific heat at constant pressure. This yields 
dT 

pCv—=$ v +pe mic -V-{F + F I ). (B7) 

To relate pressure, density and temperature, we need an equation of state for the fluid, the simplest of which is given by the 
ideal gas law, 

"~ -pT , (B8) 



where fee is Boltzmann's constant, p the mean molecular weight and mu the mass of the hydrogen atom. 

The energy generation rate per unit mass, and the opacity are often represented as power laws of temperature and density 
(Uryu & Eriguchi 1995), i.e. 

enuc = eop 7 T 5 , (B9) 

k = K p a T-t 3 , (BIO) 

where eo, 7, S, Ko, a and j3 are constants. For simplified models for lower main-sequence stars, appropriate values for these 
constants are 

7~ 1 , <5~4.5 , (Bll) 
and for Kramers' opacity law, 

a = l, (3 = 3.5 . (B12) 

Inclusion of the thermodynamic equation introduces temperature as another variable into the problem along with other 
variables which depend on p, p, u and T. This leads altogether to six equations and six unknowns. It also requires an additional 
surface boundary condition, which can be determined from an atmosphere calculation where the effects of external irradiation 
(if present) may also be included, just as in the case for standard stellar-structure calculations (e.g. Kippenhahn, Weigert & 
Hofmeister 1967). These equations have to be solved simultaneously to yield the solution at each time step. 



APPENDIX C: THE CALCULATION OF THE POTENTIAL 

In this appendix we derive analytical formulae for the derivatives of the gradient of the potential and describe the technique 
used in their implementation. 



CI The derivatives of the potential 

Differentiating equations (^) to (|3s]) with respect to r, we obtain 

! = m=-i 



dr 



(i + l)^m, x , 1 dC lm (r) ( _ lnim , . , tdD^ir) 



r l+2 



dr 



+ lr D (r)+r 



Or 



dr 



<m'Y ,m *{6\ct,')p{r',6',ci>') 



Hence, 



(CI) 

(C2) 
(C3) 
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1 dC' m (r) idD lm (r) „ 
+ r = 



r l + l Q r Q r 

and the r derivative of the potential becomes 



d$(r,9, 
dr 



I 



1=0 m=—l 
The derivative is given by 



-6 (r) + Zr £> (r) 



r i+2 



i 



" G E2]TT E 



dY lm {6 



1 T C lm {r)+r l D lm (r) 



1=0 m=—l 
Using the definition of the spherical harmonics, we find 



d$(r,8, 



' G Y.WTl E 



21 



m= ^p;""(cose) 

Using one of the recurrence relations for Legendre polynomials 

d m (cose)] 



-j--C lm (r) +r l D lm (r) 



I (cost/) = (—1) sin 
we then obtain 



d(cos6) r 



n m+lf o\ i -,\m + l ■ m+1 a d m [P(COS^)] 

P (cos6>) = (-1) sin 61 — — — 

v ' y ' d(cos0) m+1 



/ -i \ m + 1 • m+1 

= (—1) sin i 



-1 \ d 
sine) d9 



P m (cos< 



d \Pi (cos 9)] _™, „. / mcostfsin 
■ - P (cos 9) 



(-l) m sin m 6> 

m — 1 , 



dO 



Rearranging these, we get the required relation 



sin m t 



d[Pi m {cos9)} 



i (cos 



m cos 



P, m (cos0) . 



d9 y ' sm9 

When m = /, we set P m+1 equal to zero. The 9 derivative is then 



d$(r,9, 



^ 21 + 1 ^ v ' 



l 



1 = m=-l 

Finally, the (p derivative is simply given by 



priest 



cos } mcosf 



— C (r)+rD (r) 



d$(r,9, 

dcj> 



dY l 



G E^fl E ™Y lm (e,<f>) [-±-C lm (r) + r l D lm (r) 



m= — l 
1 



dm 



1 s~ilm f \ t Z j-\hn i \ 

-C (r) + r D (r) 



r i+ 



(C4) 



(C5) 



(C6) 



(C7) 



(C8) 



(C9) 



(CIO) 



(Cll) 



(C12) 



C2 Implementation 



To implement the expansion of the potential in spherical harmonics, we follow the technique described in Miiller & Stein- 
metz (1995), in which the integrals are split into a sum of integrals over sub-intervals. If we denote the position at which 
equations © & © are to be evaluated as r n , C lm {r n ) as Cl™, and D lm (r n ) as D l ™ : 



0k fBj 



j=l k=l L" ^k- 1 " 8 i- 1 

0k i-flj 



sin 9d9d<t>Y lm *(9,4>)y^ 



^' l = EE 

j = l k=l L J 0k-l J «j 



drr l+2 p(r, 9, < 



sin 9d9dcj>Y lm * (9,c/>) ^ 

i = n + l Jr i-l 



drr 1 l p(r,t 



(C13) 
(C14) 



Introducing 4™k an d By^, 



\ Lm 



0k- 1 J0 i-1 Jr i-1 



S-Ik= / / / sin 9d9d(l>drr 1 - l Y lm * (9, <j>)p(r, 9, 4>) , 

J<l>k-1 "'Sj-l Jri-i 

we may write 

c " m =EEE = ^-i + E E ^ ' 
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(C15) 
(C16) 



sin 9d9dcj)drr l+2 Y lm * (9, 0)p(r, 9, <j>) , 



=1 j=i k=i 



j=l k=l 



^=EZE = D »+* + E E • 

i=n+l j=l k=l 

We implement this as 



j=i k=i 



k )r' +2 dr , 



£2v = &4> [' sm9Y lm "{9,^)d9 [' p(r,0j 

™M = H f sm8Y lm *(8,fa)d8 /' p(r,9 h 4>^~ l dr 



(C17) 
(C18) 

(C19) 
(C20) 



For the 9 integration we use a Romberg method. A numerical integration is required due to the large 6 size of the cells near 
the poles which is a result of the equal spacing in cos 9. For the radial integration we assume that the density varies linearly 
with r within each cell 



p{r,0j,(f>k) = p(ri_i,0j,0k) + 
= Pi-ij,k + 



(r — n-i)Sp 
(n-n-i) 
(r - n_i)(5p 



(n - n-i) 



where 

5p = pi,j, k - Pi-i,j,k 
Hence 



A'" 



J,k 



Ik = <ty / 



sin^y im *(0,0 k ) 



sin 8d8Y lm * (6»,0 k ) < 



Pi-i,j,k - 



<5p 



r i+4 



n-n-i / Z + 3 (n-n-i)i + 4 





[(- 


l,j,k 






l,j,k 




[(- 


l,j,k 



l 5p 



r i _ 1 Sp 



hi r + ■ 



-3-1 



(ri— ri_i) 3— i 

2-i "r" (n-n-i) 



^Pi-i,j,k n-ri-i) 2-1 "r" (n-n-i 



^ lnr 

-3-1 



) 3-i 



V 1 



if Z = 2 , 
if Z = 3 , 
otherwise. 



(C21) 
(C22) 

(C23) 

(C24) 

(C25) 



In our calculations we consider spherical harmonics up to and including the I = 14 term. This choice follows from the results 
presented by Muller & Steinmetz (1995) who find (as we do) that the inclusion of terms of higher order does not significantly 
affect the calculation for geometries in which the use of spherical coordinates is appropriate. 



C3 Calculation of the surface derivatives 

We use the potential to calculate the surface derivatives by considering an incremental change in the potential d<E> 
9$ 9$ 9$ 

d *=8F dr+ W dB+ 8$ d +- (C26) 
It is simple to show that 

dR de ae Ir dR 



89 





SS 1 


de 


ae | R 


as 




dr 


R 







a* 1 


d(f> 








<)<!' 






dr 


R 



(C27) 



where 3> s denotes the surface potential and all terms are evaluated at the surface. When the surface is an equipotential the 
first term in the numerator is zero. 
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APPENDIX D: FINITE DIFFERENCE FORM OF EQUATIONS 

In this appendix we give the finite difference form for the terms used in equations ([Tl|) & (|l6|). As velocities are calculated on 
cell boundaries, and pressure and density at the center of cells, it is simple to represent the pressure term (from equation 
in finite difference form as 

i^p _ 2iy r (pi + ij, k -pi,j, k ) , m s 
pdr i?j,k(pi+i,j,k + Pi,j,k) 

and similarly for the derivatives in the and (j> directions. The Coriolis term in the equation of motion is solved using 
equation (flq). With the notation u,v,w to represent u r ,lt0,u</> for ease of reading, the three components of equation ( ffs| ) are 

2n< + i,j, k+ i sinflj , (D2) 



St 



St 



= 2Qw a . i , , i cos , (D3) 



i i k — w n . 



(wr_i,j, k _isin0j+^ j+ i !k _i cos0^ . (D4) 



St 

The integral form of the continuity equation (|l6|) in finite difference form is 

' ' k St ~ 5 ( ¥ ) ^~ cos0)<^ + p'+^ k (u iJj5 n-r + u i+ i J+ i k n-0 + w i+ i J k+ 1 n-0)rf + 1 J>k 5<M(- cos0)- 



i 2 ,j, k x («i-i,j,kn-r + , j+ i, k n-0 + i Jik+ i n-0)rf_ i Jjk <5^(- cos 0) + p. . + \ ^ij+i.k sin0 J+ ir i|j+ i k <5r<5(/> 

t-f- ■2— f | ^ f | ^ 

-^ i jA,k Wi 'j' kSine i-^ r i,j-ik 5r ^ + P ij! k+l Wi 'j' k + ir iJ.k+i 5r,5 (~ COs6 ') ^^i,j,k-I Wi 'j' kr iJ.k-i ,5r ' 5 (~ COs6 ') = . (D5) 



